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We theoretically and numerically study the quantum dynamics of two degenerate optical parametric oscilla¬ 
tors with mutual injections. The cavity mode in the optical coupling path between the two oscillator facets is 
explicitly considered. Stochastic equations for the oscillators and mutual injection path based on the positive P 
representation are derived. The system of two gradually pumped oscillators with out-of-phase mutual injections 
is simulated, and its quantum state is investigated. When the incoherent loss of the oscillators other than the 
mutual injections is small, the squeezed quadratic amplitudes p in the oscillators are positively correlated near 
the oscillation threshold. It indicates finite quantum correlation, estimated via Gaussian quantum discord, and 
the entanglement between the intracavity subharmonic fields. When the loss in the injection path is low, each 
oscillator around the phase transition point forms macroscopic superposition even under a small pump noise. It 
suggests that the squeezed field stored in the low-loss injection path weakens the decoherence in the oscillators. 

PACS numbers: 42.50.Ar, 42.65.Yj, 42.65.Lm, 42.50.Lc 


I. INTRODUCTION 

Optics has contributed to the fundamental test and devel¬ 
opment of quantum mechanics from its dawn. It is mainly 
because photons can have a large energy scale so that they 
are robust to thermal noise and interact with matters via co¬ 
herent multi-photon processes. Especially, an optical para¬ 
metric oscillator (OPO) Qi is a good tool to investigate 
quantum-mechanical effects, since the photon pair generated 
in the phase-matched quadratic process iH is correlated due 
to the energy and momentum conservation. There has been 
vast theoretical and ex per imental literature on squeezing |0- 
M , entanglement lITlI - l^ and quantum teleportation ll23l4^ 
based on OPOs. And now, applications of this device to state- 
of-the-art technologies such as frequency combs co¬ 

herent feedback control ll^ and quantum information pro¬ 
cessing ll^ is a new and important direction. 

Inspired by the quantum simulators ||3Q| - [^ and the con¬ 
cept of the bosonic system with artificial feedback ll^ to sim¬ 
ulate the Ising spin system, we have proposed and examined 
the systems called “coherent Ising machines” which 

are networks based on optical oscillators with mutual injec¬ 
tions. In these machines, each Ising spin is emulated by the 
field of each oscillator. Also, the Ising Hamiltonian is mapped 
to the effective photonic loss of the whole system which is 
equivalent to the sum of the gain coefficients of all the oscilla¬ 
tors. When the machine is gradually pumped from below the 
oscillation threshold, it is expected that the system oscillates 
with the minimum gain balancing the loss, and the resulting 
state corresponds to the ground state of the simulated Hamil¬ 
tonian. 

We first ha ve prop osed the machine based on injection- 
locked lasers ll34 1^. However, lasers usually do not have 
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a fixed phase for their coherent fields to a certain reference, 
thus the polarization or phase states of them cannot be binary 
as the Ising spins. However, in the Ising machine based on de¬ 
generate OPOs (DOPOs) ll^ . each DOPO above the thresh¬ 
old takes one of the binary phase states, thus they can be uti¬ 
lized to simulate the Ising spins. Simulation results show that 
it can find the ground states of all the instances of the anti¬ 
ferromagnetic Ising model (the MAX-CUT problem) on cubic 
graphs with up to twenty spins. In an experimental demon¬ 
stration a network of four OPOs successfully solved a 
problem of frustrated spins. 

An open question is if such coherent Ising machines have 
quantum features, and exploit them for their computation. In 
the quantum computation framework, theoretical evidences 
promise computational powers beyond the rich of conven¬ 
tional digital computers 0381 However, the previous theo¬ 
retical models of the coherent Ising machines, introduced for 
benchmarking, are based on the semi-classical approach. Al¬ 
though such a model can deal with up to the quadratic squeez¬ 
ing process, a doubled-variable phase space representation is 
needed to consider further non-classical effects such as macro¬ 
scopic superposition and entanglement . 

In this paper, we study a fully quantum mechanical model 
for the system of two DOPOs with mutual injections using 
the positive P representation j^ . We extend the model for 
a single DOPO 11^ by adding the signal cavity mode in the 
mutual injection path between the facets of the two DOPOs. 
We recognize that the previous work on para metric oscilla¬ 
tors 1113, [l9|] and harmonic oscillators ll43l l44l] assumes local 
evanescent-type couplings. In our work, however, we con¬ 
sider a nonlocal and dissipative coupling via the full-quantum 
model, in which the auxiliary system operators, as well as the 
primary oscillator operators, are treated by the rigorous de¬ 
scription. For the simulation, a series of stochastic differential 
equations for the signal modes are derived. Additionally, it 
is verified that the adiabatic elimination of the injection path 
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results in the standard linear injection terms for the field vari¬ 
ables. Numerical simulations on the gradually pumped system 
with out-of-phase mutual injections are performed. When the 
dissipation in the mutual injection path is moderate, the intra¬ 
cavity fields of the two DOPOs near the oscillation threshold 
can be entangled due to the quantum correlation between the 
squeezed quadrature amplitudes p = (a - a^)l(2i), which is 
estimated via Gaussian quantum discord . When the 

whole system is sufficiently closed, the injection path stores a 
squeezed vacuum. In this case, the entanglement is destroyed 
by enhanced fluctuations in the anti-squeezed quadrature am¬ 
plitude X = (a a'^) 12. However, the weak fringes of the dis¬ 
tribution functions of p in the transition indicate that macro¬ 
scopic superposition of zero-phase and 7r-phase exists in the 
DOPOs even with a small nonlinear pump noise. This sug¬ 
gests that the squeezed vacuum stored in the mutual injec¬ 
tion path protects the macroscopic superposition from deco¬ 
herence iSlH]. 

This paper is organized as follows. In Sec. [Ill we de¬ 
scribe the Ito stochastic differential equations (SDEs) for the 
field variables in the system and relate our model to the semi- 
classical model previously studied. In Sec. jllll we give the 
simulation setting and review some ingredients for the simu¬ 
lated quantities and properties. In Secs. [IV]and[Vl we discuss 
the simulation result and simulation schemes for this system. 
Sec. [Vl] concludes the paper. 
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FIG. 1. (Color online) Schematic of the system, (a) The system 
comprises two DOPOs and a mutual injection path between them as a 
cavity. The two dichroic mirrors in the injection path are assumed to 
pass the pump field completely and highly reflect the signal field, (b) 
Beamsplitter interactions between the DOPO fields and the injection- 
path mode. The spatial phase of the injection-path field should be 
considered. 


Qxp(ikcZ) = cxp(-ikcZ) = -I. 


B. System Hamiltonian 


II. THEORETICAL MODEL 


A. System overview 


The Hamiltonian for the system is an extension of that for 
a single DOPO ll^ . and can be written as 

= ^free + ^int + ^pump + ^res + ^BS, ( 1 ) 


Here, we describe the system considered in this study. 
Fig. [T] (a) shows a schematic illustration of the system. It 
is composed of two DOPOs and a mutual injection path be¬ 
tween them as a cavity. The two angled mirrors in the in¬ 
jection path are assumed to be dichroic. It means that they 
can highly reflect and confine the signal field of a frequency 
ojs while entirely transmit the driving field with a frequency 
(jOd - oJp = 2(0s, where cOp is the frequency of the pump 
mode. The coherent and real driving field Sp enters each 
DOPO to excite the pump mode. It is assumed to be classi¬ 
cal and used as the phase reference. The bosonic annihilation 
and creation operators for the pump and signal modes in the 
DOPOs are denoted as (apj, al.) and (asj, a].), where j = 1,2 
is the index for the DOPOs. Also, those for the signal mode 
in the injection path are written as al). Fig. [T] (b) dis¬ 
plays the coupling between the DOPO signal fields and the 
injection field described as beamsplitter interactions. Here, 
the field in the injection path interacts with the two DOPO 
fields at distant points, thus we have to consider the spatial 
phase explicitly. The phase factors for the bosonic operators 
at the facet of DOPO#2 depend on the injection path length. 
When we set the z axis as shown in[T](a), they are written as 
cLc exp (ikcZ) and al exp (-ikcZ), where kc is the wave number 
for the signal mode. The mutual injections are in-phase cou¬ 
plings for QxpiikcZ) = exp(-/kcz) = 1, and out-of-phase for 


where the free Hamiltonian for the relevant modes is 

2 

*T~^free — ^ F fiCOgalaQ 

The quadratic nonlinear interaction Hamiltonian is 
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( 2 ) 


Bin, = F ? 

y=i 


(3) 


where k denotes the coupling coefficient in terms of the sub¬ 
harmonic mode. The driving Hamiltonian is described by 

2 

Bpump = ih ^ exp {-iojdt) - SpCipj exp (iojdt)] , (4) 


where the classical pumping flux Sp is set to be positive and 
real. The reservoir Hamiltonian for the signal and pump 
modes is written as 


^res ^ + i^Rsj^lj + ^pj^Rpj + ^Rpj^^jj) 


7=1 




( 5 ) 
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where f rsj, f Rpj and f rc are the heat bath operators for the sig¬ 
nal, pump and signal mode in the injection path. Finally, the 
beamsplitter interaction Hamiltonian between the injection- 
path mode and DOPO signal modes is denoted by 

'J~(bs = {^c^h ~ ^l^si + , ( 6 ) 

where f is the interaction coefficient. Note that "Hbs consid¬ 
ers a part of dissipation of the signal fields from the DOPO 
cavities as coherent transmissions and reflections. 


Here, a = (asu c^pi, ctpi, O'cV and/? = iPsUPsi.Pp\.Ppi, 
Pc)^ contain ten c-number variables to describe the state. \a) = 
kpi) l^^c) and ()0*| = <^*| <^; 2 l <^;il 
are the coherent product states for the total system. The pos¬ 
itive P representation gives a positive and appropriately nor¬ 
malized distribution function for every quantum state, ax and 
Px undergo statistically independent processes in probabilistic 
simulations while they are complex conjugate in average, i.e. 
{ax) = (fix)*’ Here, X is the index for the cavity modes. 


C. Stochastic equations 

With the standard technique to treat the thermal bath m, 
we have the master equation for the density operator of the 
system. Here, we neglect the thermal detuning term. Further¬ 
more, we introduce the positive P representation for the 
five modes to expand the density operator with a well-behaved 
probability distribution function P (a,fi) 

p= f P(a,p)^-^^d^°ad^°p. (7) 

J (P \(X) 

I 


We substitute Eq. (|7]) into the master equation and use the 
operator algebra IMTII for the probability distribution descrip¬ 
tion. After switching to the rotating frame with the driving 
frequency ojd for the pump and 0 Jdl 2 for the signal mode, we 
obtain the Fokker-Planck equation (FPE) for the distribution 
P(a,fi): 


jPic.in, 


E i^s)Psj |(Tp + 
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P(a,fi), 


( 8 ) 


where 6 = kcZ, and As = oJs - <^j/2 and Ap = ojp - ojd are 
the detuning between the cavity modes and the driving field. 
The components of the last two columns in Eq. ^ come from 
the beamsplitter coupling. We may choose a more convenient 
reference for 6 when the model is expanded for many-cavity 
systems. 


OJs2 


- (js + iAP as2 + K/3s2ap2 - fo'ce"’ 

fis2 


-(js- iAs)l3s2 + Kas2l3p2 - ^I3ce~'^ 


-r 


^^p2 

1/2 

■ dW^s2{t) ■ 

Ki3p2 


dWiis2{t) 


, ( 10 ) 


With the Ito’s rule lEi which gives the correspondence be¬ 
tween the EPEs and SDEs, we reach a series of Ito SDEs for 
the c-number variables a and fi 
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C^p2 


^p {7p'^^^p)^p2 2 

Pp2 


^p ~ {yp ~ ^^p)fip2 ~ 2^52 


or/ 

1/2 

dWQ;p2{t) 

/ 0 


dWpp2{t) 


( 12 ) 


D. Adiabatic elimination of the injection path 

To see the correspondence between our model and one for 
cavity systems with coherent external injections (e.g. Ref. 
[ 3 ^, we further consider the limit where the injection-path 
mode is adiabatically eliminated, i.e. ^ 7s- The field in 
the injection path at the steady state is given by 


ac 



7^. 


-{yc-i^s')Pc-^Ps\+i^Ps2e-'\ 


or/ 

1/2 

dW^cP) 

O 

I-H 


dWpcit) 


(13) 


where dWxif) is the real Wiener increment statistically inde¬ 
pendent of each other. This corresponds to the noise term 
in the equivalent Langevin equation whose autocorrelation is 
a delta function. Adding oscillators and injection paths is 
straightforward, thus this model will give a fully quantum me¬ 
chanical treatment of mutually injecting oscillator networks. 

Here, we consider the case of a resonant driving As = Ap = 
0 and zero temperature T^ = T^ = T^ = 0. In addition, we 
adiabatically eliminate the pump variables with an assumption 
that the pump fields decay sufficiently faster than the signal 
fields. We can take diagonal diffusion amplitude matrices for 
the DOPO signal fields, then we have a simplified model as 
follows 


dasi = 


K / K y \ 

- ysO^sl + — - -0^,1 psl + 


7p 


dt 


dPsi = 


yp 


d(Xs2 — 


Ts^sl 
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ri)' 
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2/^32 )^s 2 

-^PcC 


dt 


dt 


-iO 


dt 


(17) 


al^ = y[-^a,,+e'0^a,2), (20) 

= yyipsi+e-'^msi). (21) 

At the stable steady state under the phase locking due to the 
mutual injections, «a^i), </3^i)) = {{(^s 2 )e^^,{J3s2)e~^^) is ex¬ 
pected. Hence, the injection path will be empty there, i.e. 
«) = = 0. 

Substituting Eqs. (l20b and (|2T1) into ([T4l) - ([TTb . we have 
the SDEs for the intracavity signal fields with adiabatic elim¬ 
ination of the pump and the injection path. When we further 
define the effective signal loss y' and the normalized beam¬ 
splitter coupling ^ as 


r. + 


r 


Jc 


ysjc + ip ’ 


( 22 ) 


then we have the normalized SDEs for the signal modes with 
the eliminated injection path as 


dnj = [-nj + “ dj) + dT + g - rf.dW^j{T), 

(23) 

dl^j = \-l^j + Jlj - nj) + dr + g ^A-pjdW^jir). 

(24) 


Here, rjj = gUsj, jJ.] = gPsj and g = k/ is the normal- 

ized parametric gain coefficient serving as a noise parameter. 
A = Spjsth is the normalized pumping rate and Stu = 7's7plf< is 
the classical oscillation threshold. The time is scaled with the 
signal cavity lifetime, i.e. r = y'r. dWpjir) and dWp.ij) are 
rescaled Wiener increments. The linear mutual injection terms 
with a explicit coupling phase ^rfke^^ and have the 

same forms as those introduced in the semi-classical model 
|[^ . Therefore, the theoretical framework studied here guar¬ 
antees the validity of the standard injection model also in the 
quantum mechanical phase space representations, if the dy¬ 
namics in the injection path can be neglected. 


III. SIMULATION SETTING 
A. Working equations 


dac = [-7c0^c - ^c^si + ^o^s2e'^) dt, 
dPc = {-JcPc - ^PsX + ^Ps 2 e-^)dt. 


(18) Here, we discuss and review other elements important for 
the simulation. Eirst, we describe the simulation setting. In 
^ ^ this study, we focus on the out-of-phase mutual injections. 
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namely = -1, expecting the out-of-phase corre¬ 

lation between the two macroscopic DOPO fields. We nor¬ 
malize Eqs. (O - (O as 


y su'd j 

+ 

J 

1 

■7?) + 


dT + gy 

jA-rjjdW,j(T), 






(25) 

~ysnPj 

-r 

1 


^nPc 

\dT + g^ 

jA-p]dW^.{T). 






(26) 

dpc = 

i-ycnVc 

-Cnm 

-(nri2)dT, 

(27) 

dpc = 

{-ycnPc 

-(nPl 

-^„P2)dT, 

(28) 


In this study, we define the quadrature amplitudes in the 
DOPOs as Xj = (aj -t- _^)/2 and pj = (aj - a\)l(2i). The 
distribution functions of them are given by the diagonal 
element of the density operator for the corresponding eigen¬ 
states 


P{Xj) = {Xj\psj\Xj) 




{xj\asj){ji]\xj) 

PMsj) d^a,jd%j 


= J" P{asj,Psj)e 


(33) 


where 





JsJc 

Jsyc + C^' 

jI 

ysyc + ^^^ 


(29) 

(30) 

(31) 


and simulate Eqs. (l25]) - hereafter. The time unit in all 
the results is the effective cavity lifetime 1 /y'. 

Here, we refer again to the fact that the beamsplitter cou¬ 
pling in this model can explicitly take into account a large 
part of dissipation from the DOPO cavities. The rest incoher¬ 
ent decay, which may phenomenologically include absorption 
and scatting in the oscillator, is considered by the conventional 
parameter y^. For high-2 oscillators, we can expect the case 
^ > js, and this is an important condition for the system to 
show non-trivial quantum effects. We set f = 1, and y^ and jc 
can be smaller in the simulations. 

The noise parameter g = kI determines the typi¬ 

cal order of the photon number inside the DOPOs above the 
oscillation threshold. Basically, we focus on the case of poten¬ 
tially larger photon numbers from the practical point of view, 
thus fix this parameter as g ~ 0.01. This gives 1/g^ ~ 10000 
photons in the DOPOs at oscillation. For different decay pa¬ 
rameters y' and y^, the nonlinearity k is changed accordingly 
to keep the value of g. 


B. Observable moments and distribution functions 


r {pj\asj){l3]\pj) 

P(Pj) — (PjlPsjlPj) — J PiO'sj^Psj) ^= 1 = 1 ^ ^ d O'sjd /3sj 

= J" P(asj,/3sj) d^asjd^/3sj. 

(34) 


Here, psj is the partial density operator for the signal field in 
DOPO# 7 , with the other states traced out. Also, we have used 
the inner products involving the eigenstates and coherent state 

{xj\asj) = “ exp|-2x2 + 2xjasj -~ 

(PMsj) = (^) ' exp - apjOsj + ^ - > 

(36) 

(filjlasj) = exp ^ + asjPsj^ ■ (37) 

It is known that an oscillation in P{pj) is the evidence for ex¬ 
istence of superposition components of coherent states Js^]. 
In Eq. ([31, we see that in the manifold called classical sub¬ 
space where agj and jSsj are real, the oscillation in P(pj) comes 
from the integration of the component exp ^-i2pj (agj - 
Here, the quantum noise causing stochastic discrepancies be¬ 
tween asj and JSgj is found to be essential for observing the 
fringes. Also, in the other phase space representations where 
effectively JS = a*, P(Pj) does not show any fringe due to the 
real exponent in the integrand. 


Drummond and Gardiner El have shown that normally or¬ 
dered moments of the single mode oscillator can be obtained 
by the expectation value of corresponding c-number products. 
The trivial extension to the two-mode case with the corn- 
mutability of bosonic operators for spatially separate modes 
indicates 




(32) 

Here, {a} = [P] = and other irrelevant 

modes are traced out. We consider and simulate the moments 
up to second order including ones with different modes, by 
unrestricted sampling Monte Carlo integration. 


C. Criterion for entanglement 

To examine the entanglement between two intracavity sig¬ 
nal fields, we adopt the criterion proposed by Duan et al 
ll^ . Here, we consider the pair of Einstein-Podolsky-Rosen 
(EPR)-type operators ifTTI] u+ = xi X 2 , V- = p\-p 2 and their 
fluctuation operators lsu+ = u+ - {u+) and Av_ = v_ - (v-). 
The quadrature amplitudes defined here satisfy the commu¬ 
tation relation ^Xj, p]^ = idjk12. Thus, the condition for the 
entanglement (inseparability) between the two DOPO signal 
fields is given by 

(Afi^)-r (Av^) < 1. (38) 
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D. Quantum discord 


Finally, we refer to quantum correlation estimated in this 
study. It is referred to as the property of a composite sys¬ 
tem which a local measurement changes the state of the whole 
system. It is a weaker but broader characteristic than entan¬ 
glement, showing that even separable states can have some 
quantum features such as lack of complete distinguishability 
due to a nonorthogonal basis, and pure quantumness of each 
partial system. 

Quantum discord is a measure of quantum correlation, 
based on two different ways to describe the mutual informa¬ 
tion of a bipartite system. Suppose we have a system AB com¬ 
posed of partial systems A and B. The mutual information 
based on the total system entropy is 

I(Pab) = S (Pa) ^S(pb)-S (pab), (39) 

where S(p) = -Tr (p log p) is the von Neumann entropy. On 
the other hand, the mutual information based on the condi¬ 
tional entropy S (A\B) varies with the measurement basis for 
B, because the local measurement can perturb the total sys¬ 
tem. To evaluate the genuine quantum correlation, the mea¬ 
surement basis which disturbs the system least is chosen. The 
conditional mutual information in quantum theory is hence 
defined as 

J^(Pab) = S (pa) - inf V piS (pA\i), (40) 

' I 

where i is the index for the components of the POVM mea¬ 
surement basis {Ilf} for B. pA\i is the posterior state of A pro¬ 
vided that the ith state is measured at B. The quantum discord 
is defined as the difference of them 


D^(Pab) - 1(Pab) - J^(Pab)‘ (41) 


A bipartite system with a finite discord surely has quantum 
correlation between its elements. A system without entangle¬ 
ment can have nonzero discord, and it has been reported that 
such a “dirty” state niay be available for a nontrivial speedup 
in certain problems 11 ^ with a quantum computing model 
called D2C7 jsS. 

In general, the optimization about the measurement basis 
in Eq. (l40b is hard. Its time complexity has been shown to be 
NP-complete for qubit systems 11^ . and been an open prob¬ 
lem for continuous-variable states. However, for the case of 
Gaussian states and local measurement limited to Gaussian 
POVMs, analytic formulae for the discord have been 

derived. Also, it has been shown that they quantify the amount 
of genuine quantum correlation for a large part of Gaussian 
states, including two-mode squeezed states, coherent states 
and the vacuum state iBsIl . In this study, we approximate the 
signal field in each oscillator as an Gaussian state and estimate 
the discord of them as this Gaussian quantum discord. Here, 
we consider the unnormalized quadrature amplitudes for the 
two modes [f] = 2 [xi,pi, V 2 ,p 2 ]. Then, a two-mode Gaus¬ 
sian state is characterized with the covariance matrix of them 


o-G 


:^{fjh + hrj) ■ 


■ {rMh) 


7m \ 
7m Pm / 


(42) 


where Pm 7m are 2 x 2 matrices. The state can also 
be equivalently featured by the quantities called symplectic 
invariants defined as 


As = detofM, Bs = det^^, Cg = dety^, Dg = deto-G- (43) 

When we write the binary entropy function as /^(X) = 
(x -r ^ log (x -r log and the quantities 

called symplectic eigenvalues as y+ = 5 (^ - “ ^Dg^, 

A = -r -r 2Cg , the Gaussian quantum discord is given by 

D^icro) = fs ( V^-/B(v-)-/B(v+)+inf fg ( Vdete) . (44) 

Here, o"o is the measurement basis for the partial system B. 
6 is the covariance matrix for the partial system A after B 
has been locally measured. The last term in Eq. (l44l) can be 
optimized analytically within the range of Gaussian POVMs 
(adding Gaussian ancilla bits, symplectic transformations and 
a homodyne detection), yielding 114^ 

inf det e = 

0-0 

2C] +(Bs-\)(Ds-As) + 2\Cs I + (Bs - l)(Ds -As) 
(-1 + Bs )^ 

if (Ds-AsBs)'^<(l+Bs)CjiAs+Ds), 


AsBs -Cj+Ds- ^C^s + (Ds -Ag)^- 2Cl (AsBg + Ds) 


2Bs 


Otherwise. 


(45) 


In addition, a simpler formula for two-mode squeezed ther¬ 
mal states (including squeezed vacua) has been also derived 
as HI] 


1 + 

A bipartite state with T)^(o"g) > 1 always has entanglement 
between its elements. On the other hand, an entangled state 
can have a value of the quantum discord smaller than 1. 

In our simulations, the smaller value in those calculated 
with Eqs. (l44b . (l45l) and (l46l) is taken for each point to achieve 
a good approximation and avoid possible numerical instabili¬ 
ties, especially in the case of small pumping rates. 


IV. SIMULATION RESULT 

Here, we show the results of the numerical simulation of 
the system with out-of-phase mutual injections. The initial 
state is fixed to the vacuum state, i.e. a = P = {S. The system 
is gradually pumped (S^, meaning that the pumping rate A is 
slowly increased in time so that the DOPOs are continuously 
driven from the below to above of the oscillation threshold. 
This helps the DOPOs hold the state with the minimum ef¬ 
fective gain and avoid being dragged into an unstable solution 
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created by the mutual injection terms, we set the linear sched¬ 
ule for the pumping as 

Aft 

m ( 47 ) 

where Af and tf are the pump and time parameters for the 
final state. Note that the state is always transient because 
the pumping rate is continuously increased. Transient effects 
get clear when the decay rate in the injection path is small. 
However, the sweeping is sufficiently slow so that the DO- 
POs keep themselves stable. For the numerical integration on 
the SDEs, we adopt a second-order weak scheme 1^ origi¬ 
nally proposed by Kloeden and Platen 1^ . with a time step 
At = 2xl0-\ 

We consider two cases here. In the first case, the signal in¬ 
tracavity decay rate is varied under the condition that the 
loss in the injection path is larger than it. Here, we hold 
y^ = 2y^, meaning that the loss of the whole system is in¬ 
creased with those. In the second case, we change only y^, 
keeping y^ small (y^ = 0.01). Eqs. (1251) - (l28]) are numerically 
integrated with fixed parameters = 200, Af = 1.5, y^ = 100, 
f = 1 and g ~ 0.01 for these two. The results for the first and 
second cases are symbolized as (a) and (b) in the following 
figures, respectively. We take 50,000 stochastic trials to com¬ 
pute the observables. For the numerical stability, we introduce 
a boundary condition ll^ in the numerical algorithm to assure 
that the trajectories do not go out of the classical sub space 
\rij\ < Vt, \pj\ < ^[A. This is the case for a slowly pumped 
system where the injection path is nearly empty, because the 
steady solution of Eqs. (l23l) and ([24l) under the expected con¬ 
dition rfj = pj, r ]2 = is in this manifold. The condition 
means that we neglect the second harmonic generation, the 
reverse process of parametric down-conversion, by the state 
out of the boundary. 


A. Mean photon number 

Eig. [2] shows the transits of the mean photon number in 
the first DOPO. Those for the second oscillator are omitted 
because they look the same. The photon number rises, as 
the pumping rate increases linearly in time. The difference 
in the effective oscillation threshold and resulting intensity 
of the curves comes from the difference in the magnitude 
of the mutual injections. Eqs. (l22l) . (|2^ and (l24l) are con¬ 
venient to evaluate the threshold. When in the assumption 
that rji = Pi = -r ]2 = -J^i and the field of the injection 
path is eliminated, the effective classical threshold is given 
by Ath = I - 1^1. Eor the cases of y^ = (0.05, 0.1, 0.5, 1, 5) 
in Eig. [2] (a), the corresponding coupling coefficients and ap¬ 
proximate thresholds are ^ = (0.995, 0.980, 0.67, 0.33, 0.020) 
and Ath = (0.0050,0.020,0.33, 0.67,0.98), respectively. Also, 
for y, = (0.1, 0.5, 1, 5, 10) in Eig. [2](b), ^ = (0.999, 0.995, 
0.990,0.952,0.909) and Ath = (0.0010,0.0050,0.0099,0.048, 
0.091). Note that A = (1.5/200)^ for (a) and (b). A sufficiently 
closed system with mutual injections has a drastically smaller 
threshold than a single DOPO. At the same time, however. 



Normalized time r 



Normalized time r 


FIG. 2. (Color online) Transient of the intracavity photon number 
in a DOPO dependent on (a) the signal loss of the system under 
jc = 2y^, (b) the decay rate in the injection path y^. The normal¬ 
ized oscillation threshold depends on the effective mutual injection 
strength, varying with the parameters. 50,000 stochastic runs for 
each curve. 


quantum noise in the system interrupts oscillation of the DO- 
POs and leads to less photons around the threshold than those 
classically expected by the extrapolation of the linear region. 
Also, the curve for y^ = 0.05 in (a) shows a relaxation os¬ 
cillation due to the injection path. Such a dynamics suggests 
a possibility that the variables escape the classical subspace. 
The line y^ = 0.1 in (b) is also the case, although the oscilla¬ 
tion is outside the range of the figure. 


B. Correlation function of quadrature amplitudes 

Fig. [3] andlU displays the second order correlation functions 
for the quadrature amplitudes v and p. They are normalized 
with the products of the standard deviations of the relevant 
amplitudes. In Fig. O the negative correlation in xi and X 2 en¬ 
hances as the pumping rate is increased and hence the photon 
number in the DOPOs rises. It means that the system gives 
the macroscopic out-of-phase order in x due to the mutual in¬ 
jections, corresponding to the anti-ferromagnetic order of the 
most fundamental Ising model 7T = 0 -^ 10 -z 2 programmed in 
the system. The relaxation oscillation is seen in the curve of 
{X 1 X 2 ) for a small loss y^ of 0.05. Despite that the mutual in¬ 
jections negatively couple a and p, the curves of (/?ip 2 ) in Fig. 
[4|show that the instantaneous amplitudes pi and p 2 correlate 
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Normalized time r 



FIG. 3. (Color online) Time dependence of the correlation functions 
for the signal quadrature amplitude x. (a) The loss in the system with 
7c = and (b) only the loss in the injection path 7 c are varied, xj 
are nagatively and macroscopically correlated along with the oscilla¬ 
tion. 50,000 stochastic runs for each curve. 


positively, {pipi) vanishes as the photon number rises and 
iPi) = iPi) = 0. Thus, it indicates that {p\P 2 ) shows quan¬ 
tum correlation induced by the quantum noise and mutual in¬ 
jections. Opposite signs of correlation for the two quadrature 
amplitudes x and p suggests the entanglement of wavefunc- 
tions El of the two DOPO fields. We also got positive {x\X 2 ) 
and negative {piP 2 ) in the case of in-phase mutual injections, 
corresponding to the artificial ferromagnetic interaction. 

In Fig. [3](b) and|4](b), we see that a highly closed system 
( 7 c = 0 . 1 ) interrupts the formation of the correlation between 
the squeezed vacua in the DOPOs. Here, the injection path 
stores the squeezed vacuum from the DOPOs because f > 
7 c. Thus, it effectively works as an additional noise source 
to the DOPOs. 7 c ~ f gives good correlation in p, and the 
degradation of {piP 2 } due to a larger 7 c is not so significant. 
Note that broader peaks with 7 c = 5 and 10 come from larger 
oscillation thresholds, as seen in Fig. [2 


C. Total fluctuation of EPR-type operators 

Fig. [5] presents the total variance of the EPR-type operators 
varying with time. Here, (Aw+) -l- {Av^} < 1 represents the en¬ 
tanglement between the DOPO signal fields. When the loss in 
the system is small, the total fluctuation sharply rises from the 



Normalized time r 



FIG. 4. (Color online) Time dependence of the correlation functions 
for the signal quadrature amplitude p. (a) The loss in the system 
with 7 c = Ijs and (b) only jc are changed, pj are positively and 
microscopically correlated before the oscillation. 50,000 stochastic 
runs for each curve. 


vacuum level in spite of a respectable correlation in p, seen 
in Fig. m This means that the fluctuation of v_ = P\ - P 2 
falls under the vacuum level of 0.5 before oscillation, while 
that of ii+ = xi X 2 gets much larger than that due to the 
anti-squeezed noise field accumulated in the injection path. 
Both larger ys and yc denote more dissipation and less mutual 
injections, thus the curves in Fig. [5] (a) do not satisfy the en¬ 
tanglement criterion except for the small region up to r ~ 30. 
The oscillation in the curve for 7 ^ = 0.1 indicates that the 
phase locking of the coherent DOPO fields is hampered by 
the resonance of the injection path. This is also supported by 
the fluctuation in the quantum discord (Fig. Oa)). 

In Fig. [5] (b), when only 7 ^ is increased the total noise 
comes to drop clearly below the bound before oscillation. 
Thus, the system has the entanglement there. This is be¬ 
cause the damping of the field in the injection path gets faster 
while the system keeps a large amount of the mutual injec¬ 
tions. Here, a good part of the output fields from the DOPOs 
coherently inject to each other. It is known that an entangled 
state cannot be produced only with local operations and clas¬ 
sical communication (FOCC) thus the result here shows 
that the mutual injections can convey quantum information. 
As expected from the time range where the system shows the 
entanglement, it solely refiects the quantum correlation in p, 
i.e. {Av^_) < 0.5. In fact, {Au\) comes down only nearly to 
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Normalized time r 




0 20 40 60 80 100 120 140 

Normalized time r 



FIG. 5. (Color online) Time evolution of the total fluctuation in the 
EPR-type operators (Aw+) + (Av?.) dependent on (a) the loss in the 
system with jc = (b) the loss in the injection path (Aul) + 

(Av^) < 1 means the entanglement between the DOPO signal flelds. 
50000 stochastic runs for each curve. 


FIG. 6. (Color online) Quantum discord when the state is approxi¬ 
mated as a bipartite Gaussian state. Squeezing in the DOPOs below 
the threshold and the mutual injections give a large discord. Coher¬ 
ent flelds above the threshold in them and a coherent communication 
lead to a flnite discord. 50000 stochastic runs for each curve. 


0.5 (vacuum level) hence the total noise level is always larger 
than 0.5. This indicates that the entanglement is not perfect in 
the sense that the system shows classical correlation in x. 


D. Quantum discord 


the DOPOs well above the oscillation threshold. The lines 
without discord for = 0.05,5 in (a) and one with = 0.1 
in (b) have {x\X 2 ) = -0.957, -0.6036 and -0.9615, respectively. 

The ideal field state of the two DOPOs well above the 
threshold with the classical out-of-phase correlation and its 
covariance matrix are given by 


Fig. 0 shows the approximate quantum discord when the 
state is considered as a Gaussian state. It basically reflects the 
quantum correlation in pj (see the correspondence between 
Fig. |4]and[6]). When pj is squeezed and has positive corre¬ 
lation with that in the other DOPO, the system holds a rela¬ 
tively large discord. Moreover, many curves converge at a fi¬ 
nite value ~ 0.02, except for the ones with y^ = 0.05,1,5 
in Fig. [6] (a) and one with y^ = 0.1 in (b). It is worth noting 
that this finite discord does not attribute to the squeezing in the 
DOPOs as previously discussed iB for the case of squeezed 
thermal states, but the mixture of coherent states with classi¬ 
cal correlation. We have found that the variances in xj and pj 
quickly verge on 1/4 after oscillation in the data here, thus the 
states there are well described as coherent states. Also, the 
curves with the finite discord give an almost perfect correla¬ 
tion of X (injection-locking) at the final state such as {x\X 2 ) = 
-0.9999 and -1.0000. On the other hand, a low-loss injection 
path and significant dissipation can cancel out the discord for 


Pci = \-0^cl)2 2{-(^cl\ \{o^cl\^::^\-(^cl)l \<^cl)2 2{o^cl\ \{-(^clV 

(48) 


O-ipcl) = 


0 


V 



0 

1 

0 

0 


0 0 
4 ^ 2+1 0 

0 1 j 


(49) 


where ad is the real and positive amplitude of the coherent 
states in the DOPOs. We have found that the Gaussian dis¬ 
cord calculated with Eq. (l49l) verges on ~ 0.02356 for 
ad ^50, which is in a good agreement with the values in the 
simulation. Eq. (1481) clearly represents a mixture of Gaussian 
states, thus the result indicates a genuine quantum correlation 
between coherent states with a mutual injection path without 
excess noise. 
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E. Distribution functions for quadrature amplitudes 

Here, we focus on a low-loss case, where the distribution 
functions for the squeezed amplitude p around the oscillation 
threshold are off from Gaussian curves. Fig. [7] and [ 8 ] displays 
instantaneous distributions for x and p at some time points 
for js = 0.1, 7 c = 0-2. In Fig. [71 the distribution for x gets 
broadened as the pumping rate increases. The dashed lines are 
Gaussian fitting curves for each time point. We see that it has 
some deviation from the fitting curve at r = 33 and 35. This 
indicates that the system is at the onset of the macroscopic 
bifurcation in x. 

As shown in Fig. [ 8 l both P{pi) and P^pi) at pumping rates 
around the oscillation threshold come to have small fringes 
at the sides of their central peaks. The fringes survive until 
the clear bifurcation in P{xi) and F(x 2 ). On the other hand, 
they vanish when 7 ^ and 7 c are comparable to or larger than 
^ = 1. Therefore, P(pi) and P(p 2 ) suggests the existence 
of the macroscopic superposition of the zero-phase state and 
TT-phase state in a sufficiently closed two-DOPO system. For¬ 
mation of the superposition in such a slow pumping schedule 
with only a small 7 c means that the quantum noise stored in 
the injection path is essential in the formation of superposition 
components here. The injection path contains a squeezed field 
for the two oscillators, which protects a macroscopic superpo¬ 
sition state from decoherence ||i3, . It is worth noting that 

the theoretical model considered here is different from that 
in the previous studies. The side peaks in P(pi) and P(p 2 ) 
are as high as those in an even cat state | - a) -r \a) with 
Of ~ 0.9 although the state has a larger photon number than 
\a\^ = 0.81. This is because such a DOPO state does not cor¬ 
respond to a pure cat state. Note that the fringe signal will 
be a bit weaker than the ^ng optical cat states made with 
judicious techniques l^l^. Also, a larger g and a faster 
pumping schedule will give a clearer fringe due to the tran¬ 
sient effect, as the single DOPO case |[^ . 

Here, we add the extra squeezing of the intracavity DOPO 
fields which supports the effect by the mutual injections. Fig. 
[9] displays the variances of pi and p 2 versus time (i.e. the 
pumping rate). When the system is below the threshold, they 
decrease with the rise in the pumping rate. Following the os¬ 
cillation of the DOPOs, they get back to the value for a co¬ 
herent state and the vacuum state (0.25). The minimum value 
~ 0.043 is smaller than that for a single intracavity DOPO 
field id] (0.125, meaning -3 dB squeezing). It suggests that 
the mutual injections enhance the squeezing in the DOPOs. 

Compared to expectation values of observables, the con¬ 
vergence of the distribution functions (Fig. [ 8 ]) to the num¬ 
ber of sampling is slower, because the sampled points have to 
cover the whole space where the distribution can have a non- 
negligible value. Thus, we have taken 200,000 runs to draw 
the curves here. Simultaneous formation of the side peaks in 
both P(pi) and P(p 2 ) is a good indicator that the accuracy is 
not bad, because the two DOPOs obey the SDEs of the same 
form. However, numerical errors still lead to obvious negative 
values in some curves. Also, one of the p distribution func¬ 
tions is fluctuated a lot at some time points, leading to a larger 
fringe visibility and negative values. 




Quadrature amplitude X 2 


FIG. 7. (Color online) Distribution functions at different time points 
for (a) xi and (b) X 2 . The dashed lines are Gaussian fitting curves 
with o- = (6.6, 9.1, 13.5, 20.0) for r = (29, 31, 33, 35). 200,000 
trajectories are used, js - 0.1, jc - 0.2 and g - 0.01. 


V. DISCUSSION 

In this section, we discuss other theoretical schemes to sim¬ 
ulate the system considered here, the validity of the simulation 
in this study, and the possible contribution of the quantum ef¬ 
fects in the system to the performance of the coherent Ising 
machines. 


A. Other theoretical schemes 

First, we refer to the difficulty in the simulation in this study 
with other theoretical schemes. Regarding a numerical anal¬ 
ysis for an open quantum system, direct integration on the 
master equation with the Fock state basis is the most stan¬ 
dard method as investigated in the previous relevant studies 
||i3,ll8|]. It treats a series of ordinary differential equations for 
the components of the density matrix for the system. Single¬ 
shot numerical integration for them gives all the information 
of the solution, thus we do not have to repeat stochastic sim¬ 
ulations or take ensemble averages over a number of samples. 
Also, it is relatively easy to get a good accuracy in numerical 
integration of an ordinary differential equation. However, the 
basis has an infinite number of eigenstates hence we have to 
truncate some of them. Here, the more photons possible in the 








11 




FIG. 8. (Color online) Distribution functions at different time points 
for (a) Pi and (b) p 2 . Zoomed curves around the side peaks are added 
for both. 200,000 trajectories are used, = 0.1, 7c = 0.2 and 

g = 0.01. 




FIG. 9. (Color online) Variances of pi (blue curve) and p 2 (red 
curve). The two curves are almost identical due to the same form 
of the SDEs for each DOPO. 200,000 trajectories are used, = 0.1, 
yc = 0.2 and g = 0.01. 


system, the more eigenstates needed. In addition, the number 
of modes crucially affects the complexity of the simulation. 
When we consider two DOPOs and the injection path with m 
eigenstates for each, the number of components of the den¬ 
sity matrix is m^. This amounts to unrealistic numbers such 
as 1000^ and 10000^ thus the simulation with the parameters 
here will be too costly. 

Solving the Fokker-Planck equation will be useful if we can 
find a potential solution. However, it supposes a system at 
the steady state thus cannot treat the transient regime, which 
is thought to be important for a DOPO to have macroscopic 
superposition components IH,!!!]. Also, when the mutual 
injection path is explicitly considered, we will not be able to 
find a potential solution. 

The linearized analysis has been used to investigate the 
output squeezing spectra for the quadrature amplitudes of a 
DOPO and the entanglement between the output fields of 
two DOPOs coupled with evanescent coupling lllA fl^. We 
can apply this method to get some information about the intra¬ 
cavity fields of the DOPOs. However, this is also limited in the 
steady state, and basically gives a result in the frequency do¬ 
main. To consider the properties in the time domain, we have 
to integrate the noise spectra over the entire frequency space. 
Here, a well-known relation between the covariance of input 
and output bosonic operators iH is not satisfied. It is because 
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the input field for each DOPO originating from the mutual in¬ 
jections is not a coherent or vacuum state but a squeezed state. 
Straightforward application of this relation results in negative 
values in the spectrum of (Av?^). Thus, reconsideration on the 
theoretical model might be required. 

B. Accuracy and limit of simulation 

Next, we discuss the stability of the simulation in this study. 
We recognize that the positive P representation can give un¬ 
reliable results in some cases il- However, the simulation 
here is considered to be relatively stable, because the dynam¬ 
ics of the variables are mostly bounded in a finite manifold. 
In the theoretical model, the pump modes are adiabatically 
eliminated. This significantly helps DOPO fields be bounded 
as Ref. [b^ pointed. It is also supported by the fact that the 
Fokker-Planck equation for a sing le DOPO with the adiabatic 
elimination has the solution 11^ which comprises Gaussian 
components decaying exponentially in the phase space. In 
addition, a real and diagonal diffusion amplitude matrix as¬ 
sists them to be real. The classical solution and the gradual 
pumping scheme also help the variables be bounded. Thus, 
the simulations here are not thought to be largely affected by 
the instability of the dynamics in the complex phase space. 

Nevertheless, the mutual injection path can enhance the nu¬ 
merical error of the simulation, especially when its loss jc is 
quite small. As mentioned, a relaxation oscillation might be a 
sign for the unreliability of the simulation. We have not seen 
large variance in the amplified quantities such as the photon 
number, the correlation, and fluctuation for v. However, those 
for the squeezed amplitude p directly reflects the feature of 
the noise in the system. Thus, it can be difficult to acquire 
a good accuracy for them. We see that the curves for them 
with js = 0.05 in the first case and jc = 0.1 in the second 
have large fluctuations, around which the simulation might be 
unreliable ll^ . Thus, we have taken js = 0.1, = 0.2 and 

200,000 samples for the distribution functions to avoid large 
numerical errors. The critical fluctuations around the oscil¬ 
lation can also be another difficulty because the range of the 
variables here is big due to a small g. However, it is worth not¬ 
ing that the noise in p is rather decreased there, and that in v 
does not diverge in principle because the nonlinear pump de¬ 
pletion is fully taken into account here. To reduce the number 
of samples by removing unexpected correlation in the noise 
terms and the instability, the introduction of the gauge terms 
ll^ might be needed. 

C. Impact on coherent Ising machine 

Finally, we comment on the possibilities that the quantum 
features in the DOPOs contribute to the performance of the co¬ 
herent Ising machine, which is a DOPO network with mutual 
injections emulating the Ising model. We have seen differ¬ 


ent quantum effects for different parameter regimes. First, the 
system shows transient macroscopic superposition around the 
oscillation threshold when it has a small loss compared to the 
oscillator mirror loss f. The coherent Ising machine assigns 
the up and down spin states to the discrete phase states |0) and 
Itt). It is expected to oscillate in a phase configuration with the 
small total gain, which corresponds to the ground state of the 
Ising Hamiltonian. For an efficient ground-state search, the 
network has to avoid being trapped in local minima, whose 
number scales exponentially with the problem size in hard in¬ 
stances. The macroscopic superposition under gradual pump¬ 
ing holds the information of all the states in the whole system 
simultaneously and can potentially help the system search for 
the minimum gain mode at the onset of the oscillation. This 
is similar to the quantum parallelism in quantum computers. 
Second, a loss of the mutual injections comparable to f leads 
to the quantum correlation in terms of the squeezed amplitude 
p. Further investigations are required to clarify its effect on 
the performance of the machine. 


VI. CONCLUSION 

We have founded and simulated a fully quantum mechani¬ 
cal model of the system of two DOPOs with the mutual injec¬ 
tions. The field between the two DOPO facets is introduced 
as a cavity mode to treat the input-output relation and the field 
confinement in the mutual injection path. The model can be 
extended easily to the case of a larger network. We have 
shown that the linear mutual injection terms in the positive 
P representation are derived ab initio in the limit where the 
dynamics of the injection path is neglected. The detailed sim¬ 
ulation results for the case of out-of-phase coupling revealed 
that the gradually pumped system could have quantum cor¬ 
relation, entanglement and macroscopic superposition com¬ 
ponents. The quantum correlation and entanglement require 
noise-free mutual injections with moderate loss in its path. On 
the other hand, additional quantum noise stored in the low- 
loss injection path is rather essential to generate macroscopic 
superposition components in the DOPOs. This means that the 
closed injection path with squeezed vacuum inputs alleviates 
decoherence in the DOPOs. Quantum effects in such a sim¬ 
ple setup may also play a role not only in a coherent Ising 
machine but also in broader contexts such as nano- and opto- 
mechanics, circuit QED and superconducting devices. 
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